x = [-1, -0.5, -0.4, -0.3, 1, 1, -1];
y = [1, 1, 0.8, 1, 1, -1, -1];

vert = [x', y'];
edge = zeros(7, 2);
edge(:, 1) = [1:7]';
edge(:, 2) = mod(edge(:, 1), 7) + 1; 

param = createParam(0.00001, 0.03, 0.0000025, 0.1);

domain2D = createDomain(vert, edge, param);

dim = size(domain2D.vx, 1);

b = assemblyBoundaryVector(domain2D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

res = 30;
scale = 0.9;

[px1, py1] = meshgrid(0:res/2, 0:res/2);
[px2, py2] = meshgrid(0:res, 0:res);
[px3, py3] = meshgrid(0:res/2, 0:res/2);

px1 = reshape(px1, (res/2+1)*(res/2+1), 1) / res * 2 * 0.49 - 0.49 / 2;
px2 = reshape(px2, (res+1)*(res+1), 1) / res * 1.29 - 1.29 / 2;
px3 = reshape(px3, (res/2+1)*(res/2+1), 1) / res * 2 * 0.2 - 0.2 / 2;
py1 = reshape(py1, (res/2+1)*(res/2+1), 1) / res * 2 * 1.9 - 1.9 / 2;
py2 = reshape(py2, (res+1)*(res+1), 1) / res * 1.9 - 1.9 / 2;
py3 = reshape(py3, (res/2+1)*(res/2+1), 1) / res * 2 * 1.79 - 1.79 / 2;

shift = 1;

px1 = px1 - 0.75;
py1 = py1;
px2 = px2 + 0.35;
py2 = py2;
px3 = px3 - 0.4;
py3 = py3 - 0.1;

rx = [px1; px2; px3];
ry = [py1; py2; py3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[A, Z] = assemblyMatrix(domain2D);

theta = A \ b;

u = evaluateInteriorPoint(rx, ry, domain2D, theta);

F = evaluateInteriorDeformationGradient([domain2D.vx; rx],[domain2D.vy; ry], domain2D, theta);

% plotStrain( domain2D, rx, ry, 0 )

% plotDeformation(domain2D, rx, ry, b(1:2:2*dim), b(2:2:2*dim), u(:,1), u(:,2));